TP4 - Clustering

library(FactoMineR)
library(ggplot2)
library(gridExtra)
library(factoextra)
library(cluster)
library(mclust)
library(fpc)
library(phyloseq)

Description des données

Les données étudiées dans ce TP sont issues d’analyses chimiques de plusieurs vins de la même région de l’Italie mais correspondant à trois cépages différents. Les analyses chimiques ont permis de mesurer les quantités de 13 constituants contenus dans chacun de ces 3 types de vins.

Question 1

Commencez par lire les données et vérifier si la lecture s’est bien passée à l’aide des commandes suivantes.

url <- url("https://raw.githubusercontent.com/cmaugis/Stat-Avancees---PFB-Toulouse/main/Data/WineTP.txt")
Wine = read.table(url)
Wine[, 1] <- as.factor(Wine$label)
dim(Wine)
[1] 178  14
head(Wine)
  label Alcohol Malic.acid  Ash Alcalinity.of.ash Magnesium Total.phenols
1     1   14.23       1.71 2.43              15.6       127          2.80
2     1   13.20       1.78 2.14              11.2       100          2.65
3     1   13.16       2.36 2.67              18.6       101          2.80
4     1   14.37       1.95 2.50              16.8       113          3.85
5     1   13.24       2.59 2.87              21.0       118          2.80
  Flavanoids Nonflavanoid.phenols Proanthocyanins Color.intensity  Hue
1       3.06                 0.28            2.29            5.64 1.04
2       2.76                 0.26            1.28            4.38 1.05
3       3.24                 0.30            2.81            5.68 1.03
4       3.49                 0.24            2.18            7.80 0.86
5       2.69                 0.39            1.82            4.32 1.04
  OD280.OD315 Proline
1        3.92    1065
2        3.40    1050
3        3.17    1185
4        3.45    1480
5        2.93     735
 [ reached 'max' / getOption("max.print") -- omitted 1 rows ]

Question 2

Enoncé

A l’aide de la fonction table() déterminez la répartition des 178 vins en 3 cépages.

# A COMPLETER

Correction

Les 178 vins se répartissent selon les 3 cépages ainsi :

table(Wine$label)

 1  2  3 
59 71 48

Question 3

Enoncé

A l’aide des commandes suivantes, faites une ACP des données. Commentez.

respca <- PCA(Wine, scale.unit = TRUE, quali.sup = 1, graph = FALSE)
plot(respca, choix = "ind", habillage = 1, label = "none")
plot(respca, choix = "varcor")

Correction

ACP des données avec couleur selon les labels :

respca <- PCA(Wine, scale.unit = TRUE, quali.sup = 1, graph = FALSE)
plot(respca, choix = "ind", habillage = 1, label = "none")

plot(respca, choix = "varcor")

Commentaire : On repère bien les 3 cépages dans le premier plan factoriel. Grâce au deuxième graphique, on peut regarder comment les variables intiales sont corrélées aux deux premières composantes principales. On oublie dans la suite cette information pour faire la classification non supervisée des données et on se comparera avec la variable cépage.

Méthodes par partitionnement

Question 1

Enoncé

A l’aide de la fonction kmeans(), déterminez une classification en \(K=3\) groupes avec la méthode des Kmeans.

ClassifKmeans3 = kmeans(...)

Vous pouvez visualiser la classification obtenue avec la commande suivante :

fviz_cluster(ClassifKmeans3, data = Wine[, -1], ellipse.type = "norm", labelsize = 8) +
    ggtitle("")

A l’aide des commandes suivantes, comparez la classification obtenue avec la répartition en cépages.

table(Wine$label, ClassifKmeans3$cluster)
adjustedRandIndex(Wine$label, ClassifKmeans3$cluster)
classError(Wine$label, ClassifKmeans3$cluster)$errorRate

Correction

Classification en \(K=3\) groupes avec la méthode des Kmeans :

ClassifKmeans3 = kmeans(Wine[, -1], 3, nstart = 10)
summary(ClassifKmeans3)
             Length Class  Mode   
cluster      178    -none- numeric
centers       39    -none- numeric
totss          1    -none- numeric
withinss       3    -none- numeric
tot.withinss   1    -none- numeric
betweenss      1    -none- numeric
size           3    -none- numeric
iter           1    -none- numeric
ifault         1    -none- numeric

Visualisation de la classification obtenue :

fviz_cluster(ClassifKmeans3, data = Wine[, -1], ellipse.type = "norm", labelsize = 8) +
    ggtitle("")

table(Wine$label, ClassifKmeans3$cluster)
   
     1  2  3
  1 13  0 46
  2 20 50  1
  3 29 19  0
adjustedRandIndex(Wine$label, ClassifKmeans3$cluster)
[1] 0.3711137
classError(Wine$label, ClassifKmeans3$cluster)$errorRate
[1] 0.2977528

Question 2

Enoncé

A l’aide de la fonction pam(), déterminez une classification en \(K=3\) groupes avec la méthode PAM.

ClassifPAM3 = pam(...)

Visualisez la classification obtenue et comparez-la avec la répartition en cépage.

# A COMPLETER

Comparez cette classification avec celle des Kmeans :

# A COMPLETER

Correction

On cherche une classification en \(K=3\) classes avec PAM :

ClassifPAM3 <- pam(Wine[, -1], 3, metric = "euclidean")

Les médoides des \(3\) classes sont :

ClassifPAM3$medoids
    Alcohol Malic.acid  Ash Alcalinity.of.ash Magnesium Total.phenols
51    13.05       1.73 2.04              12.4        92          2.72
136   12.60       2.46 2.20              18.5        94          1.62
73    13.49       1.66 2.24              24.0        87          1.88
    Flavanoids Nonflavanoid.phenols Proanthocyanins Color.intensity  Hue
51        3.27                 0.17            2.91            7.20 1.12
136       0.66                 0.63            0.94            7.10 0.73
73        1.84                 0.27            1.03            3.74 0.98
    OD280.OD315 Proline
51         2.91    1150
136        1.58     695
73         2.78     472
ClassifPAM3$id.med
[1]  51 136  73

On peut visualiser cette classification dans le premier plan factoriel :

fviz_cluster(ClassifPAM3, data = Wine[, -1], ellipse.type = "norm", labelsize = 8) +
    ggtitle("")

On compare cette classification avec les cépages et avec celle des Kmeans :

adjustedRandIndex(ClassifPAM3$clustering, Wine$label)
[1] 0.3715001
table(ClassifPAM3$clustering, Wine$label)
   
     1  2  3
  1 46  2  0
  2 13 19 30
  3  0 50 18
adjustedRandIndex(ClassifPAM3$clustering, ClassifKmeans3$cluster)
[1] 0.9663286
table(ClassifPAM3$clustering, ClassifKmeans3$cluster)
   
     1  2  3
  1  1  0 47
  2 61  1  0
  3  0 68  0

Question 3

Enoncé

Dans cette question, on souhaite sélectionner le nombre de classes pour la méthode des Kmeans. A l’aide des commandes suivantes, comparez les différents critères.

PlotCriteres <- function(Data, Kmax) {
    Iintra = NULL
    CH = NULL
    Silhou = NULL
    for (k in 2:Kmax) {
        Classif = kmeans(Data, k, nstart = 10)
        Iintra = c(Iintra, sum(Classif$withinss)/nrow(Data))
        aux = silhouette(Classif$cl, daisy(Data))
        Silhou = c(Silhou, mean(aux[, 3]))
        CH = c(CH, calinhara(Data, Classif$cluster))
    }
    dfindice = data.frame(K = 2:Kmax, CH = CH, Iintra = Iintra, Silhouette = Silhou)
    g1 = ggplot(dfindice, aes(x = K, y = Iintra)) + geom_point() + geom_line() +
        ggtitle("Iintra") + theme(plot.title = element_text(size = 9)) + ylab("")
    g2 = ggplot(dfindice, aes(x = K, y = CH)) + geom_point() + geom_line() + ggtitle("Calinski-Harabasz") +
        theme(plot.title = element_text(size = 9)) + ylab("")
    g3 = ggplot(dfindice, aes(x = K, y = Silhouette)) + geom_point() + geom_line() +
        ggtitle("Silhouette") + theme(plot.title = element_text(size = 9)) + ylab("")
    grid.arrange(g1, g2, g3, ncol = 3)
}
PlotCriteres(Data = Wine[, -1], Kmax = 10)
gskmn <- clusGap(Wine[, -1], FUN = kmeans, nstart = 10, K.max = 10, B = 60)
plot_clusgap(gskmn) + theme(plot.title = element_text(size = 9)) + xlab("K")

Correction

PlotCriteres <- function(Data, Kmax) {
    Iintra = NULL
    CH = NULL
    Silhou = NULL
    for (k in 2:Kmax) {
        Classif = kmeans(Data, k, nstart = 10)
        Iintra = c(Iintra, sum(Classif$withinss)/nrow(Data))
        aux = silhouette(Classif$cl, daisy(Data))
        Silhou = c(Silhou, mean(aux[, 3]))
        CH = c(CH, calinhara(Data, Classif$cluster))
    }

    dfindice = data.frame(K = 2:Kmax, CH = CH, Iintra = Iintra, Silhouette = Silhou)
    g1 = ggplot(dfindice, aes(x = K, y = Iintra)) + geom_point() + geom_line() +
        ggtitle("Iintra") + theme(plot.title = element_text(size = 9)) + ylab("")
    g2 = ggplot(dfindice, aes(x = K, y = CH)) + geom_point() + geom_line() + ggtitle("Calinski-Harabasz") +
        theme(plot.title = element_text(size = 9)) + ylab("")
    g3 = ggplot(dfindice, aes(x = K, y = Silhouette)) + geom_point() + geom_line() +
        ggtitle("Silhouette") + theme(plot.title = element_text(size = 9)) + ylab("")
    grid.arrange(g1, g2, g3, ncol = 3)
}
PlotCriteres(Data = Wine[, -1], Kmax = 10)

gskmn <- clusGap(Wine[, -1], FUN = kmeans, nstart = 10, K.max = 10, B = 60)
plot_clusgap(gskmn) + theme(plot.title = element_text(size = 9)) + xlab("K")

Avec l’inertie intra, on peut choisir 4 ou 7 classes. Avec Silhouette, on choisit 2 classes.

Question 4

Enoncé

Etudiez la classification retenue avec le critère Silhouette.

ClassifSilh <- kmeans(...)
aux <- silhouette(ClassifSilh$cl, daisy(Wine[, -1]))
fviz_silhouette(aux) + theme(plot.title = element_text(size = 9))
# A COMPLETER

Correction

ClassifSilh <- kmeans(Wine[, -1], 2)
aux <- silhouette(ClassifSilh$cl, daisy(Wine[, -1]))
fviz_silhouette(aux) + theme(plot.title = element_text(size = 9))
  cluster size ave.sil.width
1       1  123          0.68
2       2   55          0.60

table(ClassifSilh$cluster, Wine[, 1])
   
     1  2  3
  1  9 67 47
  2 50  4  1
adjustedRandIndex(ClassifSilh$cluster, Wine[, 1])
[1] 0.3694075
fviz_cluster(ClassifSilh, data = Wine[, -1], ellipse.type = "norm", labelsize = 8) +
    ggtitle("")

Question 5

Enoncé

En adaptant des commandes des questions précédentes, déterminez une classification des vins à partir des coordonnées de l’ACP à l’aide de la méthode des Kmeans.

DataACP = respca$ind$coord
dim(DataACP)
# A COMPLETER

Correction

DataACP = respca$ind$coord
PlotCriteres(DataACP, Kmax = 10)

gskmn <- clusGap(DataACP, FUN = kmeans, nstart = 10, K.max = 10, B = 60)
plot_clusgap(gskmn) + theme(plot.title = element_text(size = 9)) + xlab("K")

ClassifKmeansPCA = kmeans(DataACP, 3)
fviz_cluster(ClassifKmeansPCA, data = Wine[, -1], ellipse.type = "norm", labelsize = 8) +
    ggtitle("")

table(Wine[, 1], ClassifKmeansPCA$cluster)
   
     1  2  3
  1 59  0  0
  2  3 65  3
  3  0  0 48
adjustedRandIndex(Wine[, 1], ClassifKmeansPCA$cluster)
[1] 0.897495

Classification hiérarchique

Question 1

Enoncé

Comparez les classifications en 3 classes obtenues avec une CAH avec le lien de ward etle lien complet respectivement.

d <- dist(Wine[, -1], method = "euclidian")
hward <- hclust(....)
clward <- cutree(...)
hcomplete <- hclust(...)
clcomplete <- cutree(...)
adjustedRandIndex(clward, clcomplete)
table(clward, clcomplete)
fviz_dend(hward, k = 3, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
    show_labels = FALSE, labels_track_height = 0.8)
fviz_dend(hcomplete, k = 3, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
    show_labels = FALSE, labels_track_height = 0.8)

Correction

Kfixe = 3
d <- dist(Wine[, -1], method = "euclidian")
hward <- hclust(d, method = "ward.D2")
clward <- cutree(hward, Kfixe)
hcomplete <- hclust(d, method = "complete")
clcomplete <- cutree(hcomplete, Kfixe)
adjustedRandIndex(clward, clcomplete)
[1] 0.7540837
table(clward, clcomplete)
      clcomplete
clward  1  2  3
     1 43  5  0
     2  0 47 11
     3  0  0 72
fviz_dend(hward, k = 3, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
    show_labels = FALSE, labels_track_height = 0.8)

fviz_dend(hcomplete, k = 3, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
    show_labels = FALSE, labels_track_height = 0.8)

Question 2

Enoncé

On se concentre maintenant sur la classification hiérarchique par la méthode de Ward. Tracez la courbe des poids de l’arbre et sélectionnez un nombre de classes.

ggplot(data.frame(K = 1:20, height = sort(hward$height, decreasing = T)[1:20]), aes(x = K,
    y = height)) + geom_point() + geom_line() + ylab("") + ggtitle("height") + theme(plot.title = element_text(size = 9))

Visualisez la classification obtenue et comparez-la à la répartition des cépages.

# A COMPLETER

Correction

ggplot(data.frame(K = 1:20, height = sort(hward$height, decreasing = T)[1:20]), aes(x = K,
    y = height)) + geom_point() + geom_line() + ylab("") + ggtitle("height") + theme(plot.title = element_text(size = 9))

clward4 <- cutree(hward, 4)
table(Wine[, 1], clward4)
   clward4
     1  2  3  4
  1 26 20 13  0
  2  2  0 18 51
  3  0  0 27 21
adjustedRandIndex(Wine[, 1], clward4)
[1] 0.2818747
fviz_dend(hward, k = 4, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
    show_labels = FALSE, labels_track_height = 0.8)

Question 3

Enoncé

Mettez en place une classification hiérarchique à partir des coordonnées de l’ACP. Comparez à la répartition par cépages.

# A COMPLETER

Correction

d <- dist(DataACP, method = "euclidian")
hwardacp <- hclust(d, method = "ward.D2")
ggplot(data.frame(K = 1:20, height = sort(hwardacp$height, decreasing = T)[1:20]),
    aes(x = K, y = height)) + geom_point() + geom_line() + ylab("") + ggtitle("height") +
    theme(plot.title = element_text(size = 9))

clwardacp <- cutree(hwardacp, 3)
fviz_dend(hwardacp, k = 3, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
    show_labels = FALSE, labels_track_height = 0.8)

table(Wine[, 1], clwardacp)
   clwardacp
     1  2  3
  1 59  0  0
  2  6  8 57
  3  0 47  1
adjustedRandIndex(Wine[, 1], clwardacp)
[1] 0.7605122

Modèles de mélanges

On va maintenant s’intéresser à l’utilisation de modèles de mélanges gaussiens pour déterminer une classification des vins.

Question 1

Enoncé

A l’aide des fonctions mclustBIC() et Mclust(), déterminez une classification des données en considérant

  • toutes les formes de mélanges gaussiens
  • un nombre de classes maximal à 9 utilisant
  • le critère BIC pour la sélection de modèle
# A COMPLETER

Correction

modBIC = mclustBIC(Wine[, -1], G = 1:9)
summary(modBIC)
Best BIC values:
             VVE,3       EVE,4       VVE,4
BIC      -6849.387 -6873.64376 -6885.51547
BIC diff     0.000   -24.25637   -36.12808

On récupère le modèle sélectionné :

mselectBIC = Mclust(Wine[, -1], x = modBIC)

On visualise les valeurs du critère BIC pour les différentes formes et les différentes valeurs de \(K\)

fviz_mclust_bic(mselectBIC) + theme(plot.title = element_text(size = 9))

Question 2

Enoncé

A l’aide des commandes suivantes, étudiez les probabilités conditionnelles d’appartenance.

df = data.frame(Cluster = as.factor(apply(mselectBIC$z, 1, which.max)), probapost = apply(mselectBIC$z,
    1, max))
ggplot(df, aes(x = Cluster, y = probapost)) + geom_boxplot()

Correction

df = data.frame(Cluster = as.factor(apply(mselectBIC$z, 1, which.max)), probapost = apply(mselectBIC$z,
    1, max))
ggplot(df, aes(x = Cluster, y = probapost)) + geom_boxplot()

Les vins ont tous une probabilité conditionnelle d’appartenance maximale très forte.

Question 3

Enoncé

Visualisez la classification obtenue à l’aide de la fonction fviz_mclust()et comparez-la avec la répartition par cépages.

# A COMPLETER

Correction

fviz_mclust(mselectBIC)

table(mselectBIC$classification, Wine[, 1])
   
     1  2  3
  1 59  0  0
  2  0 69  0
  3  0  2 48
adjustedRandIndex(mselectBIC$classification, Wine[, 1])
[1] 0.9667411

Question 4

Enoncé

Reprenez les questions précédentes pour la sélection de modèles avec le critère ICL.

# A COMPLETER

Correction

modICL = mclustICL(Wine[, -1], G = 1:9)
summary(modICL)
Best ICL values:
             VVE,3       EVE,4       EVE,3
ICL      -6850.694 -6879.45651 -6890.05211
ICL diff     0.000   -28.76282   -39.35842
plot(modICL)

On séelctionne donc le même modèle avec le critère ICL que le critère BIC.